1. LECTURA DE DADES

En aquest document farem el mateix estudi que amb les dades de Mallorca però utilitzant les de Menorca. Els càlculs i procediments són anàlegs al de l’illa de Mallorca, per això obviarem alguns detalls.

En primer lloc, grafiquem en forma de sèrie temporal aquestes dades:

turisme <- read_excel("IBESTAT.xls")
turisme$Data <- gsub("M","-",turisme$Data)
gastos.ts<-ts(turisme[-1], start=c(2015,10), frequency = 12)

men <- data.frame(x=1:86, y=turisme$Men) #dades de Menorca

serie_men <- ts(men$y,start=c(2015,10),frequency = 12)

plot.ts(serie_men, main="Menorca", xlab="Any", ylab="Despeses mensuals en €")

Les dades van des d’octubre de 2015 a novembre de 2022 (llavors tenim 7 cicles complets). Podem observar que l’efecte de la COVID és veu clarament en 2020.

Per veure millor l’estacionalitat visualitzem cada període mensualment:

seasonplot(serie_men, col=c("brown", "blue","red", "orange", "pink", "purple", "yellow","green"),year.labels=TRUE, main="Estacionalitat de Menorca", xlab="Mes", ylab=" Despeses en €")

legend("bottomright",
       legend = c(2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022),
       fill =  c("brown", "blue","red", "orange", "pink", "purple", "yellow","green"),      
       border = "black")

Observam que efectivament, així com es comenta en el treball, el període de confinament va ser de principi de març, on comença a decréixer el valor dels preus, fins a finals de juny, on pareix que s’han tornat a recuperar els valors anteriors a la COVID.

Dividirem la sèrie en tres trams:

  • Des de 10-2015 fins 12-2018 (serie1_men), amb el que predirem un període. És a dir predirem de 1-2019 a 12-2019. (amb aquest comprovam que els mètodes funcionen)

  • Des de 10-2015 fins 12-2019 (serie2_men), amb la que predirem el període de la COVID (és a dir, l’any 2020)

  • Des de 10-2015 fins 12-2020 (serie3_men), amb la que predirem el períodes després de la COVID (és a dir, l’any 2021)

Visualitzem-los:

serie1_men <- ts(men$y[1:39],start=c(2015,10),frequency = 12)
serie2_men <- ts(men$y[1:51],start=c(2015,10),frequency = 12)
serie3_men <- ts(men$y[1:63],start=c(2015,10),frequency = 12)

par(mfrow=c(2,2), mar=c(4,4,4,1)+.1)
plot.ts(serie1_men, main="Serie abans de la COVID \nmenys un cicle", xlab="Any", ylab="Euros")
plot.ts(serie2_men, main="Serie abans de la COVID", xlab="Any", ylab="Euros")
plot.ts(serie3_men, main="Serie abans i durant de la COVID", xlab="Any", ylab="Euros")
plot.ts(serie_men, main="Serie completa", xlab="Any", ylab="Euros")


2. AJUST DELS MODELS:

Abans d’aplicar algun model, estudiem un poc la sèrie:

serie_men.lm <- lm(y~x, data=men) 

plot.ts(men$y, main = "Menorca", ylab = "Despeses mensuals en €", xlab="Índex de cada mes")

#dibuixam la recta de regressió sobre els punts
abline(serie_men.lm, col='red')

Si dibuixam la recta de regressió sobre les nostres dades, tot i que aquestes estan molt disperses i no s’ajusten bé a la recta, podem observar que la tendència és una mica creixent, encara que és manté més o menys constant.

boxplot(serie_men~cycle(serie_men), xlab = "mesos", ylab = "Despeses en €", main="Boxplot de Menorca")

Podem observar també la presència d’estacionalitat, que pren els seus valors màxims a la temporada d’estiu i els seus mínims en l’hivern, fet que corresponen amb les dades turístiques a les illes.

Aplicarem diversos models per ajustar la nostra sèrie i fer-ne una predicció per llavors comparar quin és el millor.

Vegem com actua cada model:


2.1. ABANS DE LA COVID


- REGRESSIÓ SEGMENTADA

Com hem comentat anteriorment la recta de regressió no s’ajusta bé a les dades, de fet el valor del \(R^2\) és molt baix:

serie1_men.lm <- lm(y~x, data=men[1:39,]) 
summary(serie1_men.lm)$adj.r.squared
## [1] 0.03739615

Per això utilitzam el paquet segmented per ajustar a una recta de regressió a trossos.

Anem a utilitzar la funció selgmented() per veure quants de punts de canvi detecta:

punts_canvi_serie1_men <-selgmented(serie1_men.lm, Kmax = 10, type="bic")
## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 7 .. 8 .. 9 .. 10 .. 
## 
## BIC to detect no. of breakpoints:
##        0        1        2        3        4        5        6        7 
## 423.9598 430.9145 432.3182 431.8441 426.2749 397.3849 386.0818 390.6075 
##        8        9     <NA> 
## 392.2832 399.0984 400.0984 
## 
## No. of selected breakpoints: 6

Aplicam la funció segmented() que ens calcula la regressió segmentada:

serie1_men.seg <- segmented(serie1_men.lm, seg.Z=~x, psi = c(4,12,16,23,28,35))

Aquesta funció ens demana que introduïm els valors on es troben els punts de canvi, i aquesta ens dona el valor estimat. Els valors reals dels punts de canvi son:

summary(serie1_men.seg)$psi
##        Initial      Est.    St.Err
## psi1.x       4  4.206349 0.6612491
## psi2.x      12 11.801361 0.3959586
## psi3.x      16 15.390922 0.3696944
## psi4.x      23 22.697828 0.4490507
## psi5.x      28 28.566720 0.4342304
## psi6.x      35 35.595164 0.3499531

Que corresponen, aproximadament, a: 1-2016, 9-2016, 12-2016, 8-2017, 2-2018 i 9-2018.

Obtenim que el valor de \(R^2\) és bastant alt

summary(serie1_men.seg)$adj.r.squared
## [1] 0.8366384

Anem a visualitzar la regressió segmentada damunt les nostres dades

p_serie1_men <- ggplot(men[1:39,], aes(x=x, y=y)) + geom_line()+
  ylim(c(0,1300))+
  ggtitle('Regressió lineal i segmentada sobre les dades \nde Menorca abans de la COVID') +
  xlab('Índex del mes') +
  ylab('Despeses mensuals en €')

my.coef1_men <- coef(serie1_men.lm) #coeficients de la recta de regressió lineal

p_serie1_men <- p_serie1_men + geom_abline(intercept=my.coef1_men[1], slope=my.coef1_men[2], color="green") #afegim la recta de regressió lineal

my.model1_men <- data.frame(x=1:39, y=fitted(serie1_men.seg)) #model nou amb els valors ajustats de la regressió segmentada

p_serie1_men <- p_serie1_men + geom_line(data=my.model1_men, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines1_men <- serie1_men.seg$psi[,2]# guardam on estan els punts de canvi estimats

p_serie1_men <- p_serie1_men + geom_vline(xintercept=my.lines1_men, linetype="dashed")

p_serie1_men <- p_serie1_men + theme(panel.background = element_blank()) + #eliminam el fons i la quadrícula del gràfic 
  geom_vline(xintercept=0) + geom_hline(yintercept=0)

ggplotly(p_serie1_men)

Visualment també es pot observar que la recta de regressió segmentada s’ajusta millor a les nostres dades.

Anem a calcular les equacions d’aquestes rectes, sabem que les rectes tenen la forma \(y=mx+n\), on \(m\) és la pendent i \(n\) el valor de tall de l’eix y.

Hi ha una funció del paquet segmented que ens dona aquesta valors de \(n\):

intercept(serie1_men.seg)
## $x
##                Est.
## intercept1   824.68
## intercept2   139.10
## intercept3  3347.90
## intercept4 -1263.70
## intercept5  3720.40
## intercept6 -2821.90
## intercept7  8470.20

També tenim una altra funció que ens calcula les pendents:

slope(serie1_men.seg)
## $x
##            Est. St.Err. t value CI(95%).l CI(95%).u
## slope1  -87.227  39.496 -2.2085  -168.570   -5.8827
## slope2   75.759  16.690  4.5391    41.385  110.1300
## slope3 -196.140  39.496 -4.9660  -277.480 -114.7900
## slope4  103.490  16.690  6.2008    69.119  137.8700
## slope5 -116.090  21.112 -5.4990  -159.570  -72.6130
## slope6  112.920  16.690  6.7658    78.549  147.3000
## slope7 -204.310  39.496 -5.1729  -285.660 -122.9700

Aleshores, la nostra recta segmentada és:

\(\hat{y}= \left\{ \begin{array}{lcc} -87.227x + 824.68 & si & x \leq \psi_1 \\ \\ 75.759x + 139.10 & si & \psi_1 < x \leq \psi_2 \\ \\ -196.140x + 3347.90 & si & \psi_2 < x \leq \psi_3 \\ \\ 103.490x - 1263.70 & si & \psi_3 < x \leq \psi_4 \\ \\ -116.090x + 3720.40 & si & \psi_4 < x \leq \psi_5 \\ \\ 112.920x - 2821.90 & si & \psi_5 < x \leq \psi_6 \\ \\ -204.310x + 8470.20 & si & \psi_6 < x \\ \end{array} \right.\)

on \(\psi_1= 4.21\), \(\psi_2=11.8\), $_3=15.39 $, $_4=22.7 $, $_5=28.57 $ i \(\psi_6= 35.6\)

Vegem els errors:

accuracy(serie1_men.seg)
##                         ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -2.915759e-15 70.70987 57.10189 -1.192331 8.628459 0.2902722

Anem a fer la predicció de 1-2019 a 12-2019.

El darrer punt de canvi que tenim és en setembre de 2018 i sabem, seguint el patró obtingut amb les dades d’entrenament, que els següents és donaran en gener de 2019, setembre de 2019 i gener de 2020. Necessitam calcular les pendents de les rectes d’entre setembre de 2018 i gener de 2020, per poder fer la predicció i calcular l’error.

Per calcular les pendents de les noves rectes farem la mitjana de les pendents anteriors. De les pendents ja calculades en el model obviarem la primera i la darrera, ja que no són vàlides. Per tant,

  • La pendent de 9-2018 a 1-2019 i de 9-2019 a 1-2020 serà : (-67.495 - 78.019)/2 = -156.115

  • La pendent de 1-2019 a 9-2019 serà : (61.003 + 69.522 + 57.015)/3 = 97.39

Ara, seguint el mateix procediment que en el cas de Mallorca, els nous punts d’intersecció són: 6755.746, -3384.454 i 8783.786.

Llavors l’equació per a la predicció és:

\(\hat{y}= \left\{ \begin{array}{lcc} -156.115x + 6755.746 & si & x \leq \psi_7 \\ \\ 97.39x - 3384.454 & si & \psi_7 < x \leq \psi_8 \\ \\ -156.115x + 8783.786 & si & \psi_8 < x \\ \end{array} \right.\)

on \(\psi_7 = 40\) i \(\psi_8 = 47\).

Visualitzem-ho:

p1_serie_men <- ggplot(men, aes(x=x, y=y)) + geom_line()+
  ggtitle('Predicció de Menorca amb el model de regressió segmentada \nabans de la COVID') + 
  xlab('Índex del mes') + 
  ylab('Despeses mensuals en €')+
  ylim(c(0,1300))

my.model1_men <- data.frame(x=1:39, y=fitted(serie1_men.seg)) #model nou amb els valors ajustats de la regressió segmentada

p1_serie_men <- p1_serie_men + geom_line(data=my.model1_men, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines1_men <- serie1_men.seg$psi[,2]# guardam on estan els punts de canvi estimats

p1_serie_men <- p1_serie_men + geom_vline(xintercept=my.lines1_men, linetype="dashed") 

p1_serie_men <- p1_serie_men + geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic

p1_serie_men <- p1_serie_men + geom_abline(intercept = 6755.746, slope=-156.115, colour="green") +
                geom_abline(intercept = -3384.454, slope=97.39, colour="blue") +
                geom_abline(intercept = 8783.786, slope=-156.115, colour="orange")


ggplotly(p1_serie_men)

Calculem l’error de la predicció:

o1_men<-c(serie_men[40:51]) #observacions reals de la sèrie  de gener de 2019 a desembre de 2019

v1_men=c(40:48)
f1_men <- sapply(v1_men, function(x) 97.39*x-3384.454) #predicció de gener 2019 a setembre 2019

v2_men=c(49:51)
f2_men <- sapply(v2_men, function(x) -156.115*x + 8783.786 )#predicció de setembre 2019 a desembre 2019

p1_men<- c(f1_men,f2_men) #predicció  de 2019

sqrt(sum((p1_men-o1_men)^2)/12) #error de la predicció
## [1] 203.5378


- DECOMPOSE

Recordem la sèrie

plot.ts(serie1_men, main= 'Menorca abans de la COVID', ylab='Despeses mensuals en €',xlab='Any')

Ja hem dit anteriorment que té una tendència mínima i podem observar també que no hi ha variabilitat.

El mètode clàssic de descomposició separa la sèrie en subseries corresponents a la tendència, la estacionalitat i el renou.

Aquestes components es poden combinar de manera additiva o multiplicativa. En el nostre cas utilitzam el model additiu: \(y_t = \mu_t + S_t + a_t\)

decompose_s1_men<-decompose(serie1_men)
plot(decompose_s1_men, xlab="Any")

El decompose, per calcular aquestes noves tendències, el que fa és agafar les sis tendències anteriors i les sis posteriors de la sèrie original i en fa la mitjana. És per això que la primera que obtenim és en abril de 2016 i la darrera en juny de 2018. Notem que per a la predicció ens quedarem amb el valor de la darrera tendència del decompose.

t1_men <- decompose_s1_men$trend #tendències de la sèrie sense el tros a predir abans de la COVID
t1_men
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2015                                                                        
## 2016       NA       NA       NA 720.7354 723.6046 715.7458 707.0200 710.3858
## 2017 715.7379 724.1704 726.8771 726.2737 734.8071 744.5700 749.6183 748.9800
## 2018 753.2433 755.7308 759.9658 769.7688 777.2221 776.3413       NA       NA
##           Sep      Oct      Nov      Dec
## 2015                NA       NA       NA
## 2016 706.4712 696.1796 698.4267 706.9596
## 2017 748.7708 749.4954 748.3637 750.5696
## 2018       NA       NA       NA       NA

Els valors de les components estacionals els calcula fent la mitjana per mesos, és a dir, per calcular la componen de gener, agafa els valors de tots els geners que tenim a la sèrie original i en fa la mitjana. Per tant, només tenim 12 valors, un per a cada mes.

s1_men <- decompose_s1_men$seasonal
s1_men <-  s1_men[4:15] #estacionalitat de gener a desembre
s1_men
##  [1] -244.4097 -256.5797 -203.0405 -102.6617  115.6563  163.6352  243.0868
##  [8]  248.9930  223.4849  188.0684 -121.3293 -254.9037

Anem a fer la predicció d’aquest model

a1_men <- c(s1_men[7:12],s1_men) #estacionalitat de juliol 2018 a desembre 2019 (es per poder fer la predicció)

pred1_decompose_men <- sapply(a1_men, function(x) 776.3413 + x) #predicció de juliol  de 2018 a desembre 2019
p1_dec_men<-c(serie1_men[1:33], pred1_decompose_men)

A1_men<- data.frame("x" = men[1:51,]$x, "y" = men[1:51,]$y, "p"= p1_dec_men)

#ho dibuixam

grafica1_men_dec <- ggplot(A1_men)+
  ylim(c(300,1100))+
  geom_line(aes(x,p), color="red")+
  geom_line(aes(x,y))+
  geom_vline(xintercept=0) + geom_hline(yintercept=300)+ 
  labs(title="Predicció de Menorca abans de la COVID amb el model de descomposició", x="Índex del mes", y="Despeses mensuals en €")+ 
  theme(panel.background = element_blank())

grafica1_men_dec

Podem observar que la predicció pareix bastant bona. Calculem l’error que es comet.

sqrt(sum((c(serie_men[40:51]- pred1_decompose_men[7:18]))^2)/12) #error predicció de gener a desembre de 2019
## [1] 65.76087


Notem que també tenim una altra instrucció a R per fer prediccions d’una sèrie, la funció predict(). Aquesta és basa en un model ETS, anem a veure la predicció:

prediccio1_men <- predict(serie1_men,h=12)
plot(prediccio1_men, xlab="Any", ylab="Despeses mensuals en €")

Pareix que la predicció és bastant bona, ja que el cicle predit segueix un mateix patró que els anteriors. Anem a veure la predicció juntament amb la sèrie original:

df_prediccio1_men <- data.frame(prediccio1_men)

p1_ets_men <- data.frame("x"= 1:51, "PointForecast"=c(serie1_men,df_prediccio1_men$Point.Forecast), "Lo80" = c(rep(NA,39),df_prediccio1_men$Lo.80), "Hi80" = c(rep(NA,39),df_prediccio1_men$Hi.80), "Lo95" = c(rep(NA,39),df_prediccio1_men$Lo.95),"Hi95" = c(rep(NA,39),df_prediccio1_men$Hi.95))

grafica_pred1_ets <- ggplot((men[1:51,]))+
  geom_ribbon(data = p1_ets_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data =p1_ets_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = p1_ets_men, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció del model ETS(A,N,A) a Menorca \nabans de la COVID", x="Índex del mes", y="Despeses mensuals en €") + 
  geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank())

grafica_pred1_ets

Podem observar que la predicció és bastant bona, ja que continua seguint un mateix patró. I l’error comés és d’uns 73 euros.

accuracy(prediccio1_men,serie2_men)
##                      ME     RMSE      MAE       MPE      MAPE      MASE
## Training set   7.725199 66.13662 54.76033 -0.033064  8.005419 0.7094055
## Test set     -22.966948 72.64049 61.40002 -5.857499 10.400781 0.7954210
##                    ACF1 Theil's U
## Training set 0.05953999        NA
## Test set     0.05881637 0.5508193


- SARIMA

Anem a veure quin model obtenim considerant un model estacional.

Pel que hem vist anteriorment, podem considerar que no hi ha tendència, llavors no fa falta fer cap diferència a la part regular, no obstant, sí que cal fer una diferenciació d’orde 12. Recordem que l’ACF i el PACF son:

par(mfrow=c(1,2))
acf(serie1_men)
pacf(serie1_men)

Per a la part regular obtenim un ARIMA(1,0,2). Feim una diferenciació a la part estacional, és a dir, d’ordre 12:

serie1_men_dif <- diff(serie1_men,12)
plot(serie1_men_dif, main="Sèrie sense estacionalitat", xlab="Any", ylab="Sèrie diferenciada")

Aquesta és la sèrie sense estacionalitat ni tendència, vegem com es modifiquen l’ACF i el PACF.

par(mfrow=c(2,1))
acf(serie1_men_dif, lag=36)
pacf(serie1_men_dif,lag=36)

Per a la part estacional obtenim que P=0, D=1 i Q=0.

És a dir, obtenim un ARIMA(1,0,2)(0,1,0)[12], vegem quin model ens proposa R:

auto.arima(serie1_men)
## Series: serie1_men 
## ARIMA(0,0,0)(0,1,0)[12] 
## 
## sigma^2 = 8644:  log likelihood = -160.68
## AIC=323.37   AICc=323.53   BIC=324.66

R ens proposa un ARIMA(0,0,0)(0,1,0)[12]. Per tant les propostes de model ARIMA són:

model1_men<-arima(serie1_men, order=c(1,0,2), seasonal = list(order=c(0,1,0), period=12))

model2_men <- arima(serie1_men, order=c(0,0,0), seasonal = list( order=c(0,1,0), period=12))

Amb uns errors de 76 i 77 euros cada un.

accuracy(model1_men)
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 16.36758 76.16479 51.76796 0.7966955 7.986194 0.4518997
##                     ACF1
## Training set -0.05174254
accuracy(model2_men)
##                    ME    RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 17.76993 77.3574 53.66121 0.9457849 8.259738 0.4684264 0.07328248

Anem a fer la predicció de 2019 d’aquests models i visualitzem-les.

fc_m1_men <-forecast(model1_men, h=12)
fc_m2_men <- forecast(model2_men,h=12)
# model 1

pre1_men <- data.frame("Point Forecast" = serie1_men, "Lo 80" = rep(NA,39), "Hi 80"=  rep(NA,39), "Lo 95" =  rep(NA,39), "Hi 95" =  rep(NA,39)) #dades abans de la predicció amb NA als intervals ja que només els volem per la predicció
 
pred_m1_men <-data.frame(fc_m1_men) 

grafica_m1_men <- data.frame("x" = 1:51, "PointForecast" = c(pre1_men$Point.Forecast,pred_m1_men$Point.Forecast), "Lo80" = c(pre1_men$Lo.80, pred_m1_men$Lo.80), "Hi80"= c(pre1_men$Hi.80, pred_m1_men$Hi.80), "Lo95" = c(pre1_men$Lo.95, pred_m1_men$Lo.95), "Hi95" = c(pre1_men$Hi.95, pred_m1_men$Hi.95))

grafica1_men <- ggplot(men[1:51,])+
  geom_ribbon(data = grafica_m1_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m1_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m1_men, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció de Menorca amb el model ARIMA(1,0,2)(0,1,0)[12] \nabans de la COVID", x="Índex del mes", y="Despeses mensuals en €") + 
  geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank())

grafica1_men

#model 2

pred_m2_men <-data.frame(fc_m2_men)

grafica_m2_men <- data.frame("x" = 1:51, "PointForecast" = c(pre1_men$Point.Forecast,pred_m2_men$Point.Forecast), "Lo80" = c(pre1_men$Lo.80, pred_m2_men$Lo.80), "Hi80"= c(pre1_men$Hi.80, pred_m2_men$Hi.80), "Lo95" = c(pre1_men$Lo.95, pred_m2_men$Lo.95), "Hi95" = c(pre1_men$Hi.95, pred_m2_men$Hi.95))

grafica2_men <- ggplot(men[1:51,])+
  geom_ribbon(data = grafica_m2_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m2_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m2_men, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció de Menorca amb el model ARIMA (0,0,0)(0,1,0)[12] \nabans de la COVID", x="Índex del mes", y="Despeses mensuals en €") + 
  geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank())

grafica2_men

Vegem i estudiem els errors de la predicció:

accuracy(fc_m1_men, serie_men[40:51], h=12)
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set  16.36758 76.16479 51.76796  0.7966955 7.986194 0.4518997
## Test set     -18.96009 80.00195 54.61251 -3.5716422 9.554005 0.4767307
##                     ACF1
## Training set -0.05174254
## Test set              NA
accuracy(fc_m2_men,serie_men[40:51], h=12)
##                     ME    RMSE      MAE        MPE     MAPE      MASE
## Training set  17.76993 77.3574 53.66121  0.9457849 8.259738 0.4684264
## Test set     -19.24750 80.4762 55.42417 -3.6350339 9.727838 0.4838159
##                    ACF1
## Training set 0.07328248
## Test set             NA
checkresiduals(fc_m1_men)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,2)(0,1,0)[12]
## Q* = 5.4139, df = 5, p-value = 0.3675
## 
## Model df: 3.   Total lags used: 8
checkresiduals(fc_m2_men)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,0)(0,1,0)[12]
## Q* = 6.1129, df = 8, p-value = 0.6346
## 
## Model df: 0.   Total lags used: 8
par(mfrow=c(1,2))
qqPlot(fc_m1_men$residuals, id=FALSE, mean=mean(fc_m1_men$residuals), sd=sd(fc_m1_men$residuals))
qqPlot(fc_m2_men$residuals, id=FALSE, mean=mean(fc_m2_men$residuals), sd=sd(fc_m2_men$residuals))

shapiro.test(fc_m1_men$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m1_men$residuals
## W = 0.91595, p-value = 0.006534
shapiro.test(fc_m2_men$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m2_men$residuals
## W = 0.93545, p-value = 0.02697
lillie.test(fc_m1_men$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m1_men$residuals
## D = 0.21097, p-value = 0.0001394
lillie.test(fc_m2_men$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m2_men$residuals
## D = 0.17916, p-value = 0.002857


Resum dels errors que cometen cada un dels models en aquest cas:

reg. segmentada descomposició clàssica ETS(A,N,A) ARIMA (1,0,2)(0,1,0)[12] ARIMA (0,0,0)(0,1,0)[12]
error model 70.71 66.14 76.16 77.36
error predicció 203.76 65.76 72.64 80 80.48


2.2. DURANT LA COVID


- REGRESSIÓ SEGMENTADA

El valor d’\(R^2\) per a la regressió lineal és:

serie2_men.lm <- lm(y~x, data= men[1:51,])
summary(serie2_men.lm)$adj.r.squared
## [1] 0.01782343

Aleshores utilitzem el paquet segmented per ajustar les nostres dades a una regressió lineal segmentada i millorar aquest valor.

Vegem els punts de canvi:

punts_canvi_serie2_men <-selgmented(serie2_men.lm, Kmax = 10, type="bic")
## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 7 .. 8 .. 9 .. 10 .. 
## 
## BIC to detect no. of breakpoints:
##        0        1        2        3        4        5        6        7 
## 557.4237 564.4217 568.4880 572.6311 580.3316 571.0819 573.0862 564.6621 
##        8        9       10 
## 565.2093 524.6382 530.1464 
## 
## No. of selected breakpoints: 8  (1 breakpoint(s) removed due to small slope diff)
serie2_men.seg <- segmented(serie2_men.lm, seg.Z=~x, psi = c(5,10,14,21,29,35,40,47))
summary(serie2_men.seg)$psi
##        Initial      Est.    St.Err
## psi1.x       5  3.954453 0.9076704
## psi2.x      10 11.839275 0.4112265
## psi3.x      14 15.390922 0.3938302
## psi4.x      21 22.697828 0.4783673
## psi5.x      29 28.566720 0.4625795
## psi6.x      35 35.431934 0.3948826
## psi7.x      40 40.285772 0.3812246
## psi8.x      47 47.349153 0.4179509

Aquests són en: 1-2016, 9-2016, 12-2016, 8-2017, 2-2018, 8-2018, 1-2019 i 8-2019.

Ara el valor de \(R^2\) de la segmentada és

summary(serie2_men.seg)$adj.r.squared
## [1] 0.8254133

Anem a visualitzar la regressió segmentada damunt les nostres dades

p_serie2_men <- ggplot(men[1:51,], aes(x=x, y=y)) + geom_line()+
  ylim(c(0,1300)) +
  ggtitle('Regressió lineal i segmentada sobre les dades \nde Menorca durant la COVID') +
  xlab('índex del mes')+ 
  ylab('Despeses mensuals en €')

my.coef2_men <- coef(serie2_men.lm) #coeficients de la recta de regressió lineal

p_serie2_men <- p_serie2_men + geom_abline(intercept=my.coef2_men[1], slope=my.coef2_men[2], color="green") #afegim la recta de regressió lineal

my.model2_men <- data.frame(x=1:51, y=fitted(serie2_men.seg)) #model nou amb els valors ajustats de la regressió segmentada

p_serie2_men <- p_serie2_men + geom_line(data=my.model2_men, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines2_men <- serie2_men.seg$psi[,2]# guardam on estan els punts de canvi estimats

p_serie2_men <- p_serie2_men + geom_vline(xintercept=my.lines2_men, linetype="dashed")

p_serie2_men <- p_serie2_men + geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

ggplotly(p_serie2_men)

Per poder escriure la funció necessitam les pendents i els punts d’intersecció amb l’eix OY, que ens ho donen les següents funcions:

#PUNTS D'INTERSECCIÓ
intercept(serie2_men.seg)
## $x
##                Est.
## intercept1   841.24
## intercept2   171.78
## intercept3  3347.90
## intercept4 -1263.70
## intercept5  3720.40
## intercept6 -2821.90
## intercept7  7248.30
## intercept8 -4337.70
## intercept9  9635.50
#PENDENTS
slope(serie2_men.seg)
## $x
##            Est. St.Err. t value CI(95%).l CI(95%).u
## slope1  -97.165  66.526 -1.4605  -232.510    38.184
## slope2   72.128  14.517  4.9684    42.593   101.660
## slope3 -196.140  42.075 -4.6616  -281.740  -110.540
## slope4  103.490  17.780  5.8208    67.320   139.670
## slope5 -116.090  22.490 -5.1620  -161.850   -70.337
## slope6  112.920  17.780  6.3512    76.750   149.100
## slope7 -171.290  29.751 -5.7573  -231.820  -110.760
## slope8  116.310  17.780  6.5415    80.133   152.480
## slope9 -178.800  42.075 -4.2496  -264.400   -93.200

L’error del model de regressió segmentada és (ens interessa el RMSE):

accuracy(serie2_men.seg)
##              ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  0 75.67986 63.29412 -1.310022 9.753939 0.3068017

Anem a fer la predicció del 2020:

Recordem punts de canvi en: 1-2016, 9-2016, 12-2016, 8-2017, 2-2018, 8-2018, 1-2019, 8-2019. Aleshores, seguint el mateix criteri que abans, tenim que les noves previsions seran en 1-2020, 8-2020 i 1-2021. De la mateixa manera que a Mallorca treim les noves pendents i punts d’intercessció i obtenim que la funció pel tros predit és:

\(\hat{y}= \left\{ \begin{array}{lcc} -161.173x + 8801.12 & si & x \leq \psi_9 \\ \\ 101.212x - 4842.9 & si & \psi_9 < x \leq \psi_{10} \\ \\ -161.173x + 10637.82 & si & \psi_{10} < x \\ \end{array} \right.\)

on \(\psi_9 = 52\) i \(\psi_8 = 59\).

Visualitzem-ho:

#ho graficam

p2_serie_men <- ggplot(men, aes(x=x, y=y)) + geom_line()+
  ggtitle('Predicció de Menorca amb el model \nde regressió segmentada durant la COVID') + 
  xlab('Índex del mes') + 
  ylab('Despeses mensuals en €') +
  ylim(c(0,1300))

my.model2_men <- data.frame(x=1:51, y=fitted(serie2_men.seg)) #model nou amb els valors ajustats de la regressió segmentada

p2_serie_men <- p2_serie_men + geom_line(data=my.model2_men, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines2_men <- serie2_men.seg$psi[,2]# guardam on estan els punts de canvi estimats

p2_serie_men <- p2_serie_men + geom_vline(xintercept=my.lines2_men, linetype="dashed") #línies verticals en els punts de canvi

p2_serie_men <- p2_serie_men + geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic

p2_serie_men <- p2_serie_men + geom_abline(intercept = 8801.12, slope=-161.173, colour="green") +
                geom_abline(intercept = -4842.9 , slope=101.212, colour="blue") +
                geom_abline(intercept = 10637.82 , slope=-161.173, colour="orange")

ggplotly(p2_serie_men)

Calculem l’error de la predicció:

o2_men<-c(serie_men[52:63]) # dades reals per fer predicció del 2020

v3_men=c(52:59)
f3_men <- sapply(v3_men, function(x) 101.212*x-4842.9) #predicció de gener 2020 a agost 2020

v4_men=c(60:63)
f4_men <- sapply(v4_men, function(x) -161.173*x + 10637.82) #predicció de setembre 2020 a desembre 2020

p2_men <- c(f3_men,f4_men) #predicció de 2020

sqrt(sum((p2_men-o2_men)^2)/12)
## [1] 354.9417


- DECOMPOSE

La descomposició de la sèrie en aquest cas és la següent

decompose_s2_men <- decompose(serie2_men)
plot(decompose_s2_men, xlab="Any")

On les components de tendència són:

t2_men <- decompose_s2_men$trend #tendències de la sèrie sense el tros a predir abans de la COVID
t2_men
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2015                                                                        
## 2016       NA       NA       NA 720.7354 723.6046 715.7458 707.0200 710.3858
## 2017 715.7379 724.1704 726.8771 726.2737 734.8071 744.5700 749.6183 748.9800
## 2018 753.2433 755.7308 759.9658 769.7688 777.2221 776.3413 772.0879 771.4596
## 2019 762.4079 762.7075 763.6933 761.6817 750.9704 749.0288       NA       NA
##           Sep      Oct      Nov      Dec
## 2015                NA       NA       NA
## 2016 706.4712 696.1796 698.4267 706.9596
## 2017 748.7708 749.4954 748.3637 750.5696
## 2018 768.9025 765.1746 763.9954 762.3817
## 2019       NA       NA       NA       NA

I les estacionals:

s2_men <- decompose_s2_men$seasonal #estacionalitat
s2_men <- s2_men[4:15] #estacionalitat de gener a desembre
s2_men
##  [1] -262.68892 -267.21878 -232.48462 -120.92743  115.23893  173.50351
##  [7]  256.53538  256.21233  228.83927  226.40427  -97.99781 -275.41614

Fem la predicció:

a2_men <- c(s2_men[7:12],s2_men) #estacionalitat de juliol 2019 a desembre 2029 (es per poder fer la predicció)

pred2_decompose_men <- sapply(a2_men, function(x) 749.0288 + x) #predicció de juliol de 2019 a desembre 2020

p2_dec_men<-c(serie2_men[1:45], pred2_decompose_men)

A2_men<- data.frame("x" = men[1:63,]$x, "y" = men[1:63,]$y, "p"= p2_dec_men)

grafica2_men_dec <- ggplot(A2_men)+
  geom_line(aes(x,p), color="red")+
  geom_line(aes(x,y))+
  labs(title="Predicció de Menorca durant la COVID amb el model de descomposició", x="Índex del mes", y="Despesese mensuals en €")+
  geom_vline(xintercept = 0)+ geom_hline(yintercept = 0)+ theme(panel.background = element_blank())

grafica2_men_dec

L’error que es comet és:

sqrt(sum((c(serie_men[52:63]- pred2_decompose_men[7:18]))^2)/12)
## [1] 346.3806


Vegem, igual que abans, la predicció amb la funció predict():

prediccio2_men <- predict(serie2_men,h=12)
plot(prediccio2_men, xlab="Any", ylab ="Despeses menusals en €")

Vegem com queda la predicció sobre la sèrie original:

df_prediccio2_men <- data.frame(prediccio2_men)
p2_ets_men <- data.frame("x"= 1:63, "PointForecast"=c(serie2_men,df_prediccio2_men$Point.Forecast), "Lo80" = c(rep(NA,51),df_prediccio2_men$Lo.80), "Hi80" = c(rep(NA,51),df_prediccio2_men$Hi.80), "Lo95" = c(rep(NA,51),df_prediccio2_men$Lo.95),"Hi95" = c(rep(NA,51),df_prediccio2_men$Hi.95))

grafica_pred2_ets <- ggplot((men[1:63,]))+
  geom_ribbon(data = p2_ets_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data =p2_ets_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = p2_ets_men, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció del model ETS(A,N,A) a Menorca durant la COVID", y="Despeses mensuals en €", x="Índex del mes") + 
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0)+ theme(panel.background = element_blank())

grafica_pred2_ets

Si calculam l’error comés, aquest és d’uns 351 euros.

accuracy(prediccio2_men, serie3_men)
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set   10.85421  64.23595  53.31916 0.5545536 7.951054 0.7563633
## Test set     -213.59590 350.65163 240.79195      -Inf      Inf 3.4157744
##                   ACF1 Theil's U
## Training set 0.1249660        NA
## Test set     0.4522522         0


- SARIMA

Quin model proposam nosaltres? Vegem l’ACF i el PACF.

par(mfrow=c(1,2))
acf(serie2_men)
pacf(serie2_men)

Per a la part regular obtenim un ARIMA(4,0,2).

Observam que hi ha estacionalitat, llavors feim una diferenciació d’ordre 12.

serie2_men_diff <- diff(serie2_men,12)
plot(serie2_men_diff, main="Sèrie sense estacionalitat", xlab="Any", ylab="Sèrie diferenciada")

Ara ja no s’observa l’estacionalitat, llavors hem de fer una diferenciació D=1. Vegem quins són els nous ACF i PACF.

par(mfrow=c(1,2))
acf(serie2_men_diff,lag=36)
pacf(serie2_men_diff,lag=36)

Obtenim que P=0 i Q=0. Per tant el model que nosaltres proposam és un ARIMA(4,0,2)(0,1,0)[12].

R proposa el següent model:

auto.arima(serie2_men)
## Series: serie2_men 
## ARIMA(0,0,0)(0,1,0)[12] 
## 
## sigma^2 = 7977:  log likelihood = -230.53
## AIC=463.06   AICc=463.17   BIC=464.73

Llavors els models que tenim són:

model3_men <- arima(serie2_men, order=c(4,0,2), seasonal = list(order=c(0,1,0), period=12))

model4_men <- arima(serie2_men, order=c(0,0,0), seasonal = list( order=c(0,1,0), period=12))

I els seus errors:

accuracy(model3_men)
##                    ME     RMSE     MAE        MPE     MAPE      MASE
## Training set 8.802988 69.36035 46.9251 0.06608693 7.264296 0.4210999
##                      ACF1
## Training set -0.007604448
accuracy(model4_men)
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 9.059943 78.10244 54.07602 -0.1320548 8.605173 0.4852714
##                    ACF1
## Training set 0.02022333

Vegem les prediccions

fc_m3_men <- forecast(model3_men, h=12)
fc_m4_men <- forecast(model4_men,h=12)
#primer model

pre2_men <- data.frame("Point Forecast" = serie2_men, "Lo 80" =  rep(NA,51), "Hi 80"= rep(NA,51), "Lo 95" = rep(NA,51), "Hi 95" = rep(NA,51))

pred_m3_men <-data.frame(fc_m3_men)

grafica_m3_men <- data.frame("x" = 1:63, "PointForecast" = c(pre2_men$Point.Forecast,pred_m3_men$Point.Forecast), "Lo80" = c(pre2_men$Lo.80, pred_m3_men$Lo.80), "Hi80"= c(pre2_men$Hi.80, pred_m3_men$Hi.80), "Lo95" = c(pre2_men$Lo.95, pred_m3_men$Lo.95), "Hi95" = c(pre2_men$Hi.95, pred_m3_men$Hi.95))

grafica3_men <- ggplot(men[1:63,])+
  geom_ribbon(data = grafica_m3_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m3_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m3_men, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Prediccó  de Menorca amb el model ARIMA (4,0,2)(0,1,0)[12] \ndurant la COVID", x="Índex del mes", y="Despeses mensuals en €")+
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica3_men

#segon model

pred_m4_men <-data.frame(fc_m4_men)

grafica_m4_men <- data.frame("x" = 1:63, "PointForecast" = c(pre2_men$Point.Forecast,pred_m4_men$Point.Forecast), "Lo80" = c(pre2_men$Lo.80, pred_m4_men$Lo.80), "Hi80"= c(pre2_men$Hi.80, pred_m4_men$Hi.80), "Lo95" = c(pre2_men$Lo.95, pred_m4_men$Lo.95), "Hi95" = c(pre2_men$Hi.95, pred_m4_men$Hi.95))

grafica4_men <- ggplot(men[1:63,])+
  geom_ribbon(data = grafica_m4_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m4_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m4_men, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Prediccó  de Menorca amb el model ARIMA (0,0,0)(0,1,0)[12] \ndurant la COVID", x="Índex del mes", y="Despeses mensuals en €")+
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica4_men

Vegem i estudiem els seus errors:

accuracy(fc_m3_men,serie_men[52:63], h=12)
##                       ME      RMSE      MAE        MPE     MAPE      MASE
## Training set    8.802988  69.36035  46.9251 0.06608693 7.264296 0.4210999
## Test set     -207.871895 347.65796 256.7831       -Inf      Inf 2.3043392
##                      ACF1
## Training set -0.007604448
## Test set               NA
accuracy(fc_m4_men,serie_men[52:63], h=12)
##                       ME      RMSE       MAE        MPE     MAPE      MASE
## Training set    9.059943  78.10244  54.07602 -0.1320548 8.605173 0.4852714
## Test set     -205.551667 350.77019 256.22500       -Inf      Inf 2.2993307
##                    ACF1
## Training set 0.02022333
## Test set             NA
checkresiduals(fc_m3_men)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(4,0,2)(0,1,0)[12]
## Q* = 1.9353, df = 4, p-value = 0.7477
## 
## Model df: 6.   Total lags used: 10
checkresiduals(fc_m4_men)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,0)(0,1,0)[12]
## Q* = 4.2132, df = 10, p-value = 0.9372
## 
## Model df: 0.   Total lags used: 10
par(mfrow=c(1,2))
qqPlot(fc_m3_men$residuals, id=FALSE, mean=mean(fc_m3_men$residuals), sd=sd(fc_m3_men$residuals))
qqPlot(fc_m4_men$residuals, id=FALSE, mean=mean(fc_m3_men$residuals), sd=sd(fc_m3_men$residuals))

shapiro.test(fc_m3_men$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m3_men$residuals
## W = 0.94617, p-value = 0.02186
shapiro.test(fc_m4_men$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m4_men$residuals
## W = 0.9436, p-value = 0.01717
lillie.test(fc_m3_men$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m3_men$residuals
## D = 0.16187, p-value = 0.001878
lillie.test(fc_m4_men$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m4_men$residuals
## D = 0.17622, p-value = 0.0004126


Resum dels errors que comet cada un del models anteriors:

reg. segmentada descomposició clàssica ETS(A,N,A) ARIMA (4,0,2)(0,1,0)[12] ARIMA (0,0,0)(0,1,0)[12]
error model 75.68 64.24 69.36 78.1
error predicció 354.9417 346.38 350.65 347.66 350.77


2.3. DESPRÉS DE LA COVID


- REGRESSIÓ SEGMENTADA

Anem a fer ús del paquet segmented.

Vegem els punts de canvi que, igual que a Mallorca, en aquest cas ens indica que no n’hi ha:

serie3_men.lm <- lm(y~x, data=men[1:63,]) 
punts_canvi_serie3_men <-selgmented(serie3_men.lm, Kmax = 10, type="bic")
## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 
## (search truncated at 6 breakpoints due to increasing values of BIC) 
## 
## BIC to detect no. of breakpoints:
##        0        1        2     <NA>        3        4        5 
## 700.6072 704.9516 712.8201 713.8201 714.3654 720.6737 720.9406 
## 
## No. of selected breakpoints:  0
serie3_men.seg <- segmented(serie3_men.lm, seg.Z=~x, psi = c(5,10,16,22,29,36,39,45,55,58))
summary(serie3_men.seg)$psi
##         Initial      Est.    St.Err
## psi1.x        5  4.206349 0.7961087
## psi2.x       10 11.543480 0.5612034
## psi3.x       16 16.046980 0.4894338
## psi4.x       22 22.578797 0.5367818
## psi5.x       29 28.338906 0.5479134
## psi6.x       36 36.541203 0.3634258
## psi7.x       39 39.431366 0.3513731
## psi8.x       45 46.731891 0.4875429
## psi9.x       55 55.984004 0.2442491
## psi10.x      58 58.007390 0.3267941

Ens queden els punts de canvi a: 1-2016, 9-2016, 1-2017, 8-2017, 1-2018, 10-2018, 12-2018, 8-2019, 5-2020 i 7-2020.

Ara, el valor de \(R^2\) de la segmentada és

summary(serie3_men.seg)$adj.r.squared
## [1] 0.813196

Anem a visualitzar la regressió segmentada damunt les nostres dades

p_serie3_men <- ggplot(men[1:63,], aes(x=x, y=y)) + geom_line()+
  ylim(c(0,1300))+
  ggtitle('Regressió lineal i segmentada sobre les dades de \nMenorca després de la COVID') + 
  xlab('índex del mes')+ 
  ylab('Despeses mensuals en €')

my.coef3_men <- coef(serie3_men.lm) #coeficients de la recta de regressió lineal

p_serie3_men <- p_serie3_men + geom_abline(intercept=my.coef3_men[1], slope=my.coef3_men[2], color="green") #afegim la recta de regressió lineal

my.model3_men <- data.frame(x=1:63, y=fitted(serie3_men.seg)) #model nou amb els valors ajustats de la regressió segmentada

p_serie3_men <- p_serie3_men + geom_line(data=my.model3_men, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines3_men <- serie3_men.seg$psi[,2]# guardam on estan els punts de canvi estimats

p_serie3_men <- p_serie3_men + geom_vline(xintercept=my.lines3_men, linetype="dashed") 

p_serie3_men <- p_serie3_men + geom_vline(xintercept = 0)+ geom_hline(yintercept = 0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic

ggplotly(p_serie3_men)

Per escriure la funció a trossos necessitam les pendents i els punts d’intersecció de les rectes:

#punts d'intersecció
intercept(serie3_men.seg)
## $x
##                  Est.
## intercept1     824.68
## intercept2     139.10
## intercept3    2722.00
## intercept4   -1500.40
## intercept5    3720.40
## intercept6   -2214.10
## intercept7   12566.00
## intercept8   -4119.80
## intercept9    6348.40
## intercept10 -22225.00
## intercept11   4038.00
#pendents
slope(serie3_men.seg)
## $x
##             Est. St.Err. t value CI(95%).l CI(95%).u
## slope1   -87.227  47.552 -1.8344  -183.260    8.8053
## slope2    75.759  20.094  3.7702    35.178  116.3400
## slope3  -148.000  33.624 -4.4016  -215.900  -80.0930
## slope4   115.130  25.417  4.5297    63.802  166.4600
## slope5  -116.090  25.417 -4.5675  -167.420  -64.7620
## slope6    93.318  16.407  5.6877    60.183  126.4500
## slope7  -311.160  75.186 -4.1385  -463.000 -159.3100
## slope8   112.000  20.094  5.5738    71.421  152.5800
## slope9  -112.000  13.727 -8.1594  -139.730  -84.2810
## slope10  398.380  75.186  5.2987   246.540  550.2300
## slope11  -54.371  33.624 -1.6170  -122.280   13.5340

L’error de la regressió segmentada és:

accuracy(serie3_men.seg)
##                        ME     RMSE     MAE  MPE MAPE      MASE
## Training set 1.352864e-15 85.77708 67.2011 -Inf  Inf 0.3253227

Anem a fer la predicció per a l’any 2021. Recordem que els punts de canvi es donen a 1-2016, 9-2016, 1-2017, 8-2017, 1-2018, 10-2018, 12-2018, 8-2019, 5-2020 i 7-2020. Degut a la pertorbació del COVID, sí que hi segueix havent un punt de canvi en estiu i l’altre en hivern successivament però ara no és donen al mateix mes. Per això, per predir els següents punts de canvi agafarem el mes més freqüent. Aleshores els següents punts de canvi seran en 1-2021, 8-2021 I 1-2022.

  • La pendent de 7-2020 a 1-2021 i de 8-2021 a 1-2022 serà : -171.8125
  • La pendent de 1-2021 a 8-2021 serà : 158.9174

Els nous punts d’intersecció es calculen de forma anàloga que a Mallorca. Per a les tres noves rectes obtenim que són \(n_1\)=10851.87 , \(n_2\)=−10314.85 i \(n_3\)=13166.98.

Aleshores la predicció de 2021 serà:

\(\hat{y}= \left\{ \begin{array}{lcc} -171.813x + 10851.87 & si & x \leq \psi_{11} \\ \\ 158.9174x - 10314.85 & si & \psi_{11} < x \leq \psi_{12} \\ \\ -171.8125x + 13166.98 & si & \psi_{12} < x \\ \end{array} \right.\)

on \(\psi_9 = 64\) i \(\psi_8 = 71\).

Visualitzem-la:

p3_serie_men <- ggplot(men, aes(x=x, y=y)) + geom_line()+
  ylim(c(-200,1300)) + ggtitle('Predicció de Menorca amb el model de \nregressió segmentada després de la COVID') + xlab('índex del mes')+ ylab('Despeses mensuals en €')

my.model3_men <- data.frame(x=1:63, y=fitted(serie3_men.seg)) #model nou amb els valors ajustats de la regressió segmentada

p3_serie_men <- p3_serie_men + geom_line(data=my.model3_men, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines3_men <- serie3_men.seg$psi[,2]# guardam on estan els punts de canvi estimats

p3_serie_men <- p3_serie_men + geom_vline(xintercept=my.lines3_men, linetype="dashed") 

p3_serie_men <- p3_serie_men + geom_vline(xintercept = 0)+ geom_hline(yintercept = 0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic

p3_serie_men <- p3_serie_men + geom_abline(intercept = 10851.87, slope=-171.8125, colour="green") +
                geom_abline(intercept = -10314.85 , slope=158.9174, colour="blue") +
                geom_abline(intercept = 13166.98 , slope=-171.8125, colour="orange")

ggplotly(p3_serie_men)

Vegem quin és aquest error que es comet:

o3_men<-c(serie_men[52:63]) # dades reals per fer predicció del 2021

v5_men=c(52:59)
f5_men <- sapply(v5_men, function(x) 158.9174*x-10314.85) #predicció de gener 2021 a agost 2021

v6_men=c(60:63)
f6_men <- sapply(v6_men, function(x) -171.8125*x + 13166.98) #predicció de setembre 2021 a desembre 2021

p3_men <- c(f5_men,f6_men) #predicció de 2021

sqrt(sum((p3_men-o3_men)^2)/12)
## [1] 1977.984


- DECOMPOSE

Visualitzem la sèrie descomposada:

decompose_s3_men <- decompose(serie3_men)
plot(decompose_s3_men, xlab="Any")

Les components de tendència són:

t3_men <- decompose_s3_men$trend #tendències de la sèrie sense el tros a predir abans de la COVID
t3_men
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2015                                                                        
## 2016       NA       NA       NA 720.7354 723.6046 715.7458 707.0200 710.3858
## 2017 715.7379 724.1704 726.8771 726.2737 734.8071 744.5700 749.6183 748.9800
## 2018 753.2433 755.7308 759.9658 769.7688 777.2221 776.3413 772.0879 771.4596
## 2019 762.4079 762.7075 763.6933 761.6817 750.9704 749.0288 757.8138 761.0633
## 2020 608.6500 589.2296 571.2750 548.4712 539.5446 546.7817       NA       NA
##           Sep      Oct      Nov      Dec
## 2015                NA       NA       NA
## 2016 706.4712 696.1796 698.4267 706.9596
## 2017 748.7708 749.4954 748.3637 750.5696
## 2018 768.9025 765.1746 763.9954 762.3817
## 2019 763.7283 741.4167 680.8050 632.0783
## 2020       NA       NA       NA       NA

I les d’estacionalitat:

s3_men <- decompose_s3_men$seasonal #estacionalitat
s3_men <-  s3_men[4:15] #estacionalitat de gener a desembre

Anem a fer la predicció

a3_men <- c(s3_men[7:12],s3_men) #estacionalitat de juliol 2020 a desembre 2021 (es per poder fer la predicció)

pred3_decompose_men <- sapply(a3_men, function(x) 546.7817 + x) #predicció de juliol de 2019 a desembre 2020

p3_dec_men<-c(serie3_men[1:57], pred3_decompose_men)

A3_men<- data.frame("x" = men[1:75,]$x, "y" = men[1:75,]$y, "p"= p3_dec_men)

grafica3_men_dec <- ggplot(A3_men)+
  geom_line(aes(x,p), color="red")+
  geom_line(aes(x,y))+
  labs(title="Predicció de Menorca després de la COVID amb el model de descomposició", x="Índex del mes", y="Despeses mensualns en €")+
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica3_men_dec

Calculem l’error de la predicció

sqrt(sum((c(serie_men[64:75]- pred3_decompose_men[7:18]))^2)/12)
## [1] 218.4594


Amb la funció predict() la predicció seria la següent:

prediccio3_men <- predict(serie3_men,h=12)
plot(prediccio3_men, xlab="Any", ylab="Despeses mensuals en €")

Vegem-la juntament amb les nostres dades:

df_prediccio3_men <- data.frame(prediccio3_men)
p3_ets_men <- data.frame("x"= 1:75, "PointForecast"=c(serie3_men,df_prediccio3_men$Point.Forecast), "Lo80" = c(rep(NA,63),df_prediccio3_men$Lo.80), "Hi80" = c(rep(NA,63),df_prediccio3_men$Hi.80), "Lo95" = c(rep(NA,63),df_prediccio3_men$Lo.95),"Hi95" = c(rep(NA,63),df_prediccio3_men$Hi.95))

grafica_pred3_ets <- ggplot((men[1:75,]))+
  geom_ribbon(data = p3_ets_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data =p3_ets_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = p3_ets_men, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció del model ETS(A,N,A) a Menorca després de la COVID", x="Índex del mes", y="Despeses mensuals en €")+ 
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica_pred3_ets

Obtenim un error d’uns 117 euros.

accuracy(prediccio3_men,serie_men)
##                      ME     RMSE       MAE       MPE     MAPE      MASE
## Training set   1.544179 126.6657  89.34926      -Inf      Inf 0.7824237
## Test set     -65.427782 117.3246 106.49793 -8.062414 13.97902 0.9325931
##                    ACF1 Theil's U
## Training set 0.06833943        NA
## Test set     0.46102653  1.084369


- SARIMA

Quin model proposam nosaltres? Vegem l’ACF i el PACF:

par(mfrow=c(1,2))
acf(serie3_men)
pacf(serie3_men)

Per a la part regular obtenim p=2, q=2 i d=0.

Sí hi ha indicació d’estacionalitat, llavors feim una diferenciació d’ordre 12 i tornam a calcular l’ACF i el PACF:

serie3_men_diff <- diff(serie3_men,12)
plot.ts(serie3_men_diff, main="Sèrie sense estacionalitat", ylab="Sèrie diferenciada", xlab="Any")

par(mfrow=c(1,2))
acf(serie3_men_diff)
pacf(serie3_men_diff)

Aleshores obtenim un model ARIMA(2,0,2)(1,1,1)[12].

R proposa el següent model:

auto.arima(serie3_men)
## Series: serie3_men 
## ARIMA(0,1,2)(0,1,1)[12] 
## 
## Coefficients:
##           ma1      ma2     sma1
##       -0.2359  -0.5527  -0.4631
## s.e.   0.1334   0.1258   0.2793
## 
## sigma^2 = 22849:  log likelihood = -322.38
## AIC=652.76   AICc=653.65   BIC=660.41

Llavors tenim els següents models

model5_men <- arima(serie3_men, order=c(2,0,2), seasonal = list( order=c(1,1,1), period=12))

model6_men <- arima(serie3_men, order=c(0,1,2), seasonal = list( order=c(0,1,1), period=12))

Amb uns errors de:

accuracy(model5_men)
##                     ME     RMSE      MAE  MPE MAPE      MASE        ACF1
## Training set -8.492597 122.6242 77.91071 -Inf  Inf 0.6601132 0.003622309
accuracy(model6_men)
##                     ME   RMSE      MAE  MPE MAPE    MASE         ACF1
## Training set -12.09391 130.56 82.80016 -Inf  Inf 0.70154 -0.004851061

Vegem les prediccions

fc_m5_men <- forecast(model5_men, h=12)
fc_m6_men <- forecast(model6_men,h=12)
#grafiquem-les
#primer model

pre3_men <- data.frame("Point Forecast" = serie3_men, "Lo 80" =  rep(NA,63), "Hi 80"= rep(NA,63), "Lo 95" = rep(NA,63), "Hi 95" = rep(NA,63))

pred_m5_men <-data.frame(fc_m5_men)

grafica_m5_men <- data.frame("x" = 1:75, "PointForecast" = c(pre3_men$Point.Forecast,pred_m5_men$Point.Forecast), "Lo80" = c(pre3_men$Lo.80, pred_m5_men$Lo.80), "Hi80"= c(pre3_men$Hi.80, pred_m5_men$Hi.80), "Lo95" = c(pre3_men$Lo.95, pred_m5_men$Lo.95), "Hi95" = c(pre3_men$Hi.95, pred_m5_men$Hi.95))


grafica5_men <- ggplot(men[1:75,])+
  geom_ribbon(data = grafica_m5_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m5_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m5_men, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció de Menorca amb el model ARIMA (2,0,2)(1,1,1)[12] \ndesprés de la COVID", x="Índex del mes", y="Despeses mensuals en €")+ 
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica5_men

#segon model

pred_m6_men <-data.frame(fc_m6_men)

grafica_m6_men <- data.frame("x" = 1:75, "PointForecast" = c(pre3_men$Point.Forecast,pred_m6_men$Point.Forecast), "Lo80" = c(pre3_men$Lo.80, pred_m6_men$Lo.80), "Hi80"= c(pre3_men$Hi.80, pred_m6_men$Hi.80), "Lo95" = c(pre3_men$Lo.95, pred_m6_men$Lo.95), "Hi95" = c(pre3_men$Hi.95, pred_m6_men$Hi.95))

grafica6_men <- ggplot(men[1:75,])+
  geom_ribbon(data = grafica_m6_men, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m6_men, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m6_men, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
labs(title="Predicció de Menorca amb el model ARIMA (0,1,2)(0,1,1)[12] \ndesprés de la COVID", x="Índex del mes", y="Despeses mensuals en €")+ 
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica6_men

Vegem quin és l’error de cada model i estudiem-los

accuracy(fc_m5_men,serie_men[64:75], h=12)
##                     ME     RMSE      MAE      MPE     MAPE      MASE
## Training set -8.492597 122.6242 77.91071     -Inf      Inf 0.6601132
## Test set     76.160973 137.6361 89.64163 11.16486 12.81547 0.7595056
##                     ACF1
## Training set 0.003622309
## Test set              NA
accuracy(fc_m6_men,serie_men[64:75], h=12)
##                     ME     RMSE       MAE      MPE     MAPE     MASE
## Training set -12.09391 130.5600  82.80016     -Inf      Inf 0.701540
## Test set     218.06171 257.1519 218.06171 30.86593 30.86593 1.847569
##                      ACF1
## Training set -0.004851061
## Test set               NA
checkresiduals(fc_m5_men)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,2)(1,1,1)[12]
## Q* = 9.1279, df = 7, p-value = 0.2436
## 
## Model df: 6.   Total lags used: 13
checkresiduals(fc_m6_men)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,2)(0,1,1)[12]
## Q* = 12.831, df = 10, p-value = 0.2333
## 
## Model df: 3.   Total lags used: 13
par(mfrow=c(1,2))
qqPlot(fc_m5_men$residuals, id=FALSE, mean=mean(fc_m5_men$residuals), sd=sd(fc_m5_men$residuals))
qqPlot(fc_m6_men$residuals, id=FALSE, mean=mean(fc_m6_men$residuals), sd=sd(fc_m6_men$residuals))

shapiro.test(fc_m5_men$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m5_men$residuals
## W = 0.83562, p-value = 7.234e-07
shapiro.test(fc_m6_men$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m6_men$residuals
## W = 0.86346, p-value = 4.995e-06
lillie.test(fc_m5_men$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m5_men$residuals
## D = 0.21125, p-value = 2.001e-07
lillie.test(fc_m6_men$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m6_men$residuals
## D = 0.15808, p-value = 0.0004721


Resum dels errors que comet cada un del models anteriors:

reg. segmentada descomposició clàssica ETS(A,N,A) ARIMA (2,0,2)(1,1,1)[12] ARIMA (0,1,2)(0,1,1)[12]
error model 85.78 126.67 122.62 130.56
error predicció 1977.984 218.4594 117.32 137.64 257.15